###############################
# Analysis of conjoint data 
###############################

rm(list=ls())

sink("004_conjoint_diagnostics.text")

# load packages 
library(foreign)
library(lmtest)
library(sandwich)
library(stargazer)
library(msm)
library(Hmisc)

setwd("~/Dropbox/Romanian local elections 2016/03_analysis")

# load function that does clustered SEs
vcovCluster <- function(
  model,
  cluster
)
{
  require(sandwich)
  require(lmtest)
  if(nrow(model.matrix(model))!=length(cluster)){
    stop("check your data: cluster variable has different N than model")
  }
  M <- length(unique(cluster))
  N <- length(cluster)           
  K <- model$rank   
  if(M<50){
    warning("Fewer than 50 clusters, variances may be unreliable (could try block bootstrap instead).")
  }
  dfc <- (M/(M - 1)) * ((N - 1)/(N - K))
  uj  <- apply(estfun(model), 2, function(x) tapply(x, cluster, sum));
  rcse.cov <- dfc * sandwich(model, meat = crossprod(uj)/N)
  return(rcse.cov)
}

##############################
# Prepare data
##############################

# load data
d <- read.dta("~/Dropbox/Romanian local elections 2016/01_data/conjoint_Romania.dta")
head(d)
dim(d)

d_income = d[d$income!=99,]
d_income = d_income[order(d_income$code, decreasing=F), ]
table(d_income$income)
d_income$idnum2 = rep(1:479,each=10)

##################
# Diagnostics
##################

# Carryover effects
d$fpair = as.factor(d$pair)
d1 <- lm(outcome ~ atgen + atexp + atinti + atmita + atinte + atpol + atven +  
                   fpair + 
                   atgen*fpair + atexp*fpair + atinti*fpair + atmita*fpair + atinte*fpair + atpol*fpair + atven*fpair,
                   data=d)
d1_c <- round(coeftest(d1, vcov = vcovCluster(d1, cluster = d$idnum)),4)

# Profile order effects
d$fcandidate = as.factor(d$candidate)
d4 <- lm(outcome ~ atgen + atexp + atinti + atmita + atinte + atpol + atven +  
                    fcandidate + 
                    atgen*fcandidate + atexp*fcandidate + atinti*fcandidate + atmita*fcandidate + atinte*fcandidate +     
                    atpol*fcandidate + atven*fcandidate,
                    data=d)
d4_c <- round(coeftest(d4, vcov = vcovCluster(d4, cluster = d$idnum)),4)

# Random assigment: regress gender on attributes (regress covariate on treatment)

# urban
d2 <- lm(urban ~ as.factor(atgen) + as.factor(atexp) + as.factor(atinti) + as.factor(atmita) + as.factor(atinte) + as.factor(atpol) + as.factor(atven), data=d)
d2_c <- round(coeftest(d2, vcov = vcovCluster(d2, cluster = d$idnum)),4)

# high school
d3 <- lm(highschool ~ as.factor(atgen) + as.factor(atexp) + as.factor(atinti) + as.factor(atmita) + as.factor(atinte) + as.factor(atpol) + as.factor(atven), data=d)
d3_c <- round(coeftest(d3, vcov = vcovCluster(d3, cluster = d$idnum)),4)
        
#####################
# Regression tables
#####################

# APPENDIX I
# Carryover effects
d1_c
stargazer(d1_c,title="Carryover effects",
          type = "text",
          align=TRUE, 
          omit.stat=c("LL","ser","f"), 
          covariate.labels=c("Female","Incumbent","Threat to non-supporters","100 RON","Social assistance","Investigated",
                             "Sentenced","Renovate schools","Renovate schools and roads","High income",
                             "Pair 2","Pair 3","Pair 4","Pair 5",
                             "Female*Pair 2","Female*Pair 3","Female*Pair 4","Female*Pair 5",
                             "Incumbent*Pair 2","Incumbent*Pair 3","Incumbent*Pair 4","Incumbent*Pair 5",
                             "Threat to non-supporters*Pair 2","Threat to non-supporters*Pair 3",
                             "Threat to non-supporters*Pair 4","Threat to non-supporters*Pair 5",
                             "100 RON*Pair 2","Social assistance*Pair 2","100 RON*Pair 3","Social assistance*Pair 3",
                             "100 RON*Pair 4","Social assistance*Pair 4","100 RON*Pair 5","Social assistance*Pair 5",
                             "Investigated*Pair 2","Sentenced*Pair 2","Investigated*Pair 3","Sentenced*Pair 3",
                             "Investigated*Pair 4","Sentenced*Pair 4","Investigated*Pair 5","Sentenced*Pair 5",
                             "Renovate schools*Pair 2","Renovate schools and roads*Pair 2",
                             "Renovate schools*Pair 3","Renovate schools and roads*Pair 3",
                             "Renovate schools*Pair 4","Renovate schools and roads*Pair 4",
                             "Renovate schools*Pair 5","Renovate schools and roads*Pair 5",
                             "High income*Pair 2","High income*Pair 3","High income*Pair 4","High income*Pair 5"),
          no.space=TRUE,  
          table.placement = "H",
          #single.row = TRUE,
          star.cutoffs = c(0.05,0.01,0.001),
          dep.var.caption  = "Outcome",
          dep.var.labels   = "Electoral Choice",
          out="diagnostics_carryover.htm")
                                       

# APPENDIX H
# Profile order effects
d4_c
stargazer(d4_c,title="Profile order effects",
          type = "text",
          align=TRUE, 
          omit.stat=c("LL","ser","f"), 
          covariate.labels=c("Female","Incumbent","Threat to non-supporters","100 RON","Social assistance","Investigated",
                             "Sentenced","Renovate schools","Renovate schools and roads","High income",
                             "Candidate 2",
                             "Female*Candidate 2","Incumbent*Candidate 2","Threat to non-supporters*Candidate 2",
                             "100 RON*Candidate 2","Social assistance*Candidate 2","Investigated*Candidate 2",
                             "Sentenced*Candidate 2","Renovate schools*Candidate 2","Renovate schools and roads*Candidate 2",
                             "High income*Candidate 2"),
          no.space=TRUE,  
          table.placement = "H",
          #single.row = TRUE,
          star.cutoffs = c(0.05,0.01,0.001),
          dep.var.caption  = "Outcome",
          dep.var.labels   = "Electoral Choice",
          out="diagnostics_profile.htm")


# APPENDIX F
# Random assigment:urban 
d2_c
stargazer(d2_c,title="Balance test: Urban",
          type = "text",
          align=TRUE, 
          omit.stat=c("LL","ser","f"), 
          covariate.labels=c("Female","Incumbent","Threat to non-supporters","100 RON","Social assistance","Investigated",
                             "Sentenced","Renovate schools","Renovate schools and roads","High income"),
          no.space=TRUE,  
          table.placement = "H",
          #single.row = TRUE,
          star.cutoffs = c(0.05,0.01,0.001),
          dep.var.caption  = "Outcome",
          dep.var.labels   = "Urban",
          out="diagnostics_random_urban.htm")

# APPENDIX G
# Random assigment:high school
d3_c
stargazer(d3_c,title="Balance test: High-school",
          type = "text",
          align=TRUE, 
          omit.stat=c("LL","ser","f"), 
          covariate.labels=c("Female","Incumbent","Threat to non-supporters","100 RON","Social assistance","Investigated",
                             "Sentenced","Renovate schools","Renovate schools and roads","High income"),
          no.space=TRUE,  
          table.placement = "H",
          #single.row = TRUE,
          star.cutoffs = c(0.05,0.01,0.001),
          dep.var.caption  = "Outcome",
          dep.var.labels   = "High-school",
          out="diagnostics_random_highschool.htm")

sink()
